data = readtable('../../data/main_VAR/US_full.xlsx');
date = data.year;

startdatenum = 1947;
enddatenum = 2020;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% convenience yield
load('../../data/supplement/US_convyield.mat');
cy = m.spread(2:end);

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data((data.year >= startdatenum - 1), :);
data = data((data.year <= enddatenum), :);

taxrevgdp = data.tau(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.tau(1); % government tax revenue to GDP ratio
spendgdp = data.g(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.g(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.x(2:end); % real gdp growth (in logs)
inflation = data.pi(2:end); % growth in GDP price deflator (in logs)

ynom1 = log(1 + data.x1yrCMT(2:end) / 100 + cy * 0.9); % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = log(1 + data.x10yrCMT(2:end) / 100 + cy * 0.5); % nominal yield on a 10-year Treasury note, expressed per annum. 

pdm = data.pdm(2:end);
divgrm = data.divgrm(2:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Read in debt/gdp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gdebt = data.debt_gdp_sargenthall;

% adjust for convenience yield, where m.mult accounts for how much convenience yield each maturity bucket earns
cy_gdp_a = gdebt(2:end) .* m.mult(2:end) .* m.spread(2:end);
taxrevgdp = taxrevgdp + cy_gdp_a;

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1;
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% log pd ratio on stocks
divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


tts(1) = taxrevgdp(1); 
gts(1) = spendgdp(1); 

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end

disp('mean surplus')

100 * [mean(tts - gts) mean(surplusgdp)]

run ../../tools/run_var_estimate;

save('MAT/res_annual_VAR_nodebt_US_post1946_cy.mat','-regexp', '^(?!(option|Globaloption)$).');
